#include "cstdio"

using namespace std;

namespace nglib{
    #include "nglib.h"
}

#include <meshing.hpp>

#include <cstring>

void SwapBytes(char *array, int size, int n)
{
  char *x = new char[size];
  for(int i = 0; i < n; i++) {
    char *a = &array[i * size];
    memcpy(x, a, size);
    for(int c = 0; c < size; c++) a[size - 1 - c] = x[c];
  }
  delete[] x;
}

using namespace netgen;
using namespace nglib;


//Msh your msh
//forexample ./.vtk
void WriteVTK(const Mesh* Msh,const string path){

   FILE* fp = fopen(path.c_str(),"wb");

  fprintf(fp, "# vtk DataFile Version 2.0\n");
  fprintf(fp, "%s, Created by Gmsh\n", path.c_str());
  fprintf(fp, "BINARY\n");
  fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
  fprintf(fp,"POINTS %d double\n",Msh->GetNP());


  int np = Msh->GetNP();
  int i, j, k;

  //二进制输出
  for (i = 1; i <= np; i++)
  {
      double xyz[] = {Msh->Point(i)(0),Msh->Point(i)(1),Msh->Point(i)(2)};
      //bigend
      SwapBytes(reinterpret_cast<char*>(&xyz[0]),sizeof(double)*3,3);
      fwrite(&xyz[0],sizeof(double),3,fp);   
  }
  fprintf(fp,"\n");

    int ne = Msh->GetNE();
    //计算总共的点数
    for(j=0,i=1;i<=ne;i++){
        const Element& e = Msh->VolumeElement(i);
        j+= e.GetNP();
    }

    int nse = Msh->GetNSE();
    for(i=1;i<=nse;i++){
        const Element2d& e2 = Msh->SurfaceElement(i);
        j+= e2.GetNP();
    }

//     //获取Info
//     vtk<<"CELLS "<<ne+nse<<" "<<j<<endl;
cout<<"ne+nse"<<ne+nse<<"j is"<<j+ne+nse<<endl;
    //sum of numPoints and connectivity indices over each cell
    fprintf(fp,"CELLS %d %d\n",ne+nse,j+ne+nse);//调试发现

//write 2d
    for(i=1;i<=nse;i++){
        const Element2d& e2 = Msh->SurfaceElement(i);
        int l = e2.GetNV();
        int verts[60];
        verts[0] = l;
        for(j=1;j<=l;j++){
            // vtk<<e2.PNum(j)-1<<" ";
            verts[j] = e2.PNum(j)-1;
        }
        SwapBytes(reinterpret_cast<char*>(&verts[0]),sizeof(int),l+1);
        fwrite(&verts[0],sizeof(int),1+l,fp);
    }
// //write 3d
    for(i=1;i<=ne;i++){
        const Element& e = Msh->VolumeElement(i);
        int l = e.GetNV();
        int verts[60];
        verts[0] = l;
        for(j=1;j<=l;j++){
            // vtk<<e2.PNum(j)-1<<" ";
            verts[j] = e.PNum(j)-1;
        }
        SwapBytes(reinterpret_cast<char*>(&verts[0]),sizeof(int),l+1);
        fwrite(&verts[0],sizeof(int),1+l,fp);
    }

    fprintf(fp,"\n");

    fprintf(fp,"CELL_TYPES %d\n",ne+nse);

 //write 2d
    for(i=1;i<=nse;i++){
        const Element2d& e = Msh->SurfaceElement(i);

        int type = -1;

        if(e.GetType() == SEGMENT){
            // vtk<<3<<endl;
            type = 3;
        } 
        else if(e.GetType() == TRIG){
            // vtk<<5<<endl;
            type = 5;
        }
        else if(e.GetType() == TET){
            // vtk<<10<<endl;
            type = 10;
        }else{
            //  vtk<<1<<endl;
             type = 1;
        }

        SwapBytes(reinterpret_cast<char*>(&type),sizeof(int),1);
        fwrite(&type,sizeof(int),1,fp);
    }
//write 3d
    for(i=1;i<=ne;i++){
        const Element& e = Msh->VolumeElement(i);
      
        int type = -1;

        if(e.GetType() == SEGMENT){
            // vtk<<3<<endl;
            type = 3;
        } 
        else if(e.GetType() == TRIG){
            // vtk<<5<<endl;
            type = 5;
        }
        else if(e.GetType() == TET){
            // vtk<<10<<endl;
            type = 10;
        }else{
            //  vtk<<1<<endl;
             type = 1;
        }

        SwapBytes(reinterpret_cast<char*>(&type),sizeof(int),1);
        fwrite(&type,sizeof(int),1,fp);
    }

    // fprintf(fp,"\n");
    fclose(fp);

}
